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I. INTRODUCTION 

The possibility of creating optical lattices in trapped 
Bose-condensed gases has provided an opportunity to 
study superfluids in novel situations. The presence of 
the lattice leads to a variety of solid state effects associ- 
ated with the coherent motion of the atoms in a periodic 
potential. For example, the oscillation frequency of the 
centre of mass motion of the condensate is reduced 
as a result of the enhanced effective mass of the atoms 
tunnelling between potential wells. Furthermore, when 
subjected to a uniform force, as provided by gravity or 
alternatively by accelerating the optical lattice itself 0, 
Bloch oscillations of the condensate have been observed. 
Reducing the amplitude of the lattice potential leads to 
a breakdown of these oscillations as a result of Landau- 
Zener tunnelling between bands . All of these observa- 
tions are essentially a manifestation of the superfluidity 
of the Bose condensate in an optical lattice. Another as- 
pect of equal interest is the breakdown of superfluidity as 
recently observed in a study of the centre of mass mo- 
tion of a trapped condensate moving through an optical 
lattice. When the amplitude of the oscillation exceeded 
a critical value, dissipation was seen to set in. 

In this paper we study a Bose-condensed gas subjected 
to a uniform optical lattice in a regime where the dynam- 
ics of the condensate wave function is well described by 
the time-dependent Gross-Pitaevskii (GP) equation. In 
particular, we are concerned with small amplitude col- 
lective modes which at long wavelength are phonon-like 
excitations. The relevant physical parameters determin- 
ing the properties of the excitation are the optical po- 
tential amplitude, Voj the lattice constant, d, the mean 
density, n, of the gas, and the magnitude of the super- 
current. The problem has been addressed theoretically 
in a number of papers using a variety of techniques and 
approximations. 

Berg-S0rensen and M0lmer ^ were the first to inves- 
tigate phonon excitations within an optical lattice. They 
solved the Bogoliubov equations numerically for a one- 
dimensional model and established that the long wave- 
length excitations are phonon-like, having an energy dis- 
persion that is linear in the wave vector of the mode. 
They also obtained an analytic expression for the sound 



speed, s, which is based on a combined weak potential 
and slowly varying approximation. These calculations 
were extended by Wu and Niu ^ to the case in which the 
condensate carries a current. This work is noteworthy for 
having pointed out that the modes exhibit both energetic 
and dynamic instabilities for sufficiently large currents. 
The former instability is associated with the Landau cri- 
terion for the breakdown of superfluidity, while the latter 
is related to the onset of dissipation as observed in j^. 
The recent paper by Machholm et al. explores these in- 
stabilities further. Similar results were obtained by Bron- 
ski et al. ^] by considering a special form of the lattice 
potential, while Konotop and Salerno Q used a different 
approach to establish that the dynamic (or modulational) 
instability leads to the generation of solitons. 

When the potential wells are sufficiently deep, the con- 
densate is well-localized on each lattice site and a tight- 
binding description becomes useful. Javanainen [Tcij used 
this picture within a many-body formulation to derive the 
phonon dispersion throughout the Brillouin zone for a 
one dimensional lattice. This calculation is in fact equiv- 
alent to one based on the Bogoliubov equations or the 
discrete version of the time-dependent GP equation • 
The virtue of these methods is that they provide analyt- 
ical expressions for the dispersion relation, although an 
accurate a priori determination of the tight-binding pa- 
rameters involves further numerical calculation. Smerzi 
et al. extended the results of Javanainen [T^ by de- 
riving the phonon dispersion for a current-carrying state 
and found a dynamic instability that is responsible for a 
so-called 'superfluid-insulator' transition. 

The phonon dispersion at long wavelengths was ad- 
dressed in |13| by deriving an energy functional involv- 
ing density and phase fluctuations which vary slowly in 
space. The approach is closely allied to the effective mass 
approximation used in solid state physics Il4l and re- 
cently applied to Bose gases in optical lattices [13 , and to 
multiple-scale analysis 0, 0| ■ The phonon sound speed 
is found to be 

n djjL 
V m* dn 

where m* is an effective mass and /i is the chemical po- 
tential. The precise definition of the effective mass as 
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(m*)~^ = hr'^d'^e/ dk"^ , where e is the energy per particle 
of the iV-particle system, seems first to have appeared 
later 0, ■ This is an important point since there are 
other plausible candidates for the effective mass. In fact, 
we shall see that a somewhat different, but equivalent, 
definition can be given. In this regard, we note that the 
effective mass appearing in both the effective mass jisj l 
and multiple-scale i9] theories is that corresponding to 
the bare optical potential. In other words, the effect of 
the interactions on this parameter is not included, and 
therefore the use of the dynamical equations obtained in 
these theories will not in general give the correct Bogoli- 
ubov sound speed. 

Our purpose in this paper is to obtain the long wave- 
length phonon dispersion directly from the Bogoliubov 
equations defining the collective modes. This is achieved 
by developing a systematic expansion of these equations 
in powers of the phonon wave vector q. We do this first 
for the current-free state (Sec. IIII|I . confirming the re- 
sult for the sound speed given above. We then consider a 
current-carrying state fSec. lIVII and obtain the analogous 
phonon dispersion in this case, reproducing the result ob- 
tained in I?! by means of a hydrodynamic analysis. Our 
expansion technique can be viewed as a justification of 
the assumptions on which the hydrodynamic approach is 
based. Furthermore, it provides explicit perturbative ex- 
pressions for the various physical quantities that appear 
in the theory (for example, the effective mass). 

In Sec. ^ we present the theoretical background re- 
quired for the calculation of small amplitude collective 
excitations in an optical lattice. For the most part we 
consider a three dimensional optical potential with cu- 
bic symmetry, although we also touch on systems with 
one dimensional modulation as well as radially confined 
condensates. The underlying periodicity of the the opti- 
cal potential implies that the Bogoliubov equations ad- 
mit solutions having a Bloch function form. This aspect 
accounts for the use of a Bloch function basis in solv- 
ing these equations in both the current-free (Sec. IIII|) 
and current-carrying (Sec. IIV|I states. However, differ- 
ent calculational methods are used in the two cases and 
these are therefore presented separately. We also exam- 
ine various physical limits (Thomas-Fermi, weak poten- 
tial, weak coupling and tight-binding) in order to make 
contact with previous work. As stated previously, our 
main result for the phonon dispersion affirms the result 
which follows from the insightful use of hydrodynamic 
equations to describe the dynamics of long wavelength 
fiuct nations |i7i.il3i|. 



II. BASIC THEORY 

We consider an extended 3D BEG subjected to stand- 
ing wave light fields that give rise to a periodic exter- 
nal potential having the property T4pt (r -I- R) = I^pt (r) , 
where R is a Bravais lattice vector. For the most 
part we restrict ourselves to a cubic lattice for which 



R = d(nix -t- + n^i), with m an integer. 

We base our analysis on the time-dependent Gross- 
Pitaevskii (GP) equation for the condensate wave func- 
tion, *(r,i), 

^?l|vI,(r, t)=(^-^ + yopt(r) + g|*(r, t)p) vl,(r, t). 

(1) 

This equation admits stationary solutions of the form 

*(r,t) = $(r)e-'^*/'* 
where <I>(r) satisfies the time- independent GP equation 

V2$ + Kpt$-t-g|$|2$ = ^$ (2) 



2m 

with the normalization 



d^r |$(r)|2 = iV, 



where N is the total number of particles in the volume 
V . Also of interest is the total energy of the system given 
by 



2m 



Kpt 



$ + ^ / 
2 Jv 



(3) 

The energy parameter /i is the chemical potential and is 
related to E^tot by /i = dEtot/dN. 

Often the ground state solution of the GP equation is 
of interest but we shall also consider states which have a 
superfiuid flow. These states have a Bloch function form 



$nk(r) 



where n is a band index and k is a wave vector restricted 
to the first Brillouin zone. The factor ^/H, where n is the 
mean density, is introduced in the definition of Wnk so as 
to give the normalization 



^^^c?V|u;„k(r) 



1 , 



(4) 



where f2 is the Wigner-Seitz volume. The condensate 
density is then ^^(r) |'I>„k(r)P = n|w„k(r)P- The 
Bloch function w„k(r) is in general complex and is the 
self-consistent solution of 



2 m 



V^opt + gn\Wnk\ Wnk = fJ-nkWnk 



(5) 

We assume that Wnk has the periodicity of the lattice, 
w„k(r -I- R) = u'„k(r), although it should be noted that 
period-doubled states also exist 0. The chemical po- 
tential, Unkin), is implicitly a function of the mean den- 
sity and depends on the particular Bloch state being con- 
sidered. 
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The superfluid current density in this state is 



L{r) = ^($:kV$nk-v$:k'i'nk) 

= —hk\wnk{r)f + 7^ [w*kVw„k - (Vw*k)Wnk] 
m Zmi 

and has the property V • js(r) = 0. Introducing the 
superfluid velocity according to the relation 



(6) 



we have 



and the GP equation is expanded to first order in the 
deviation (5<l>(r, t). By writing 



J$(r, t) = Mj(r)e" 



-iEit/h _ y*(j.yiE'it/ti 



(12) 



where Ei is allowed to be complex, one obtains the fol- 
lowing Bogoliubov equations for the quasiparticle ampli- 
tudes Ui and Wi, 

Lu,{y) - g^l^{Y)v,{Y) = E,u,{r) 

Lv,{v) ^ g^*X{v)u,{v) = -E,v,{v) (13) 



V3(r) = _(k + V0„k(r)) 

m 



(7) 



where 6'„k(r) is the phase of the Bloch function Wn\^. The 
spatially-averaged superfluid velocity is 



i / dVv,(r) = -(k+(V0„k)). 



(8) 



In one dimension, the periodicity of w„j,(a;) implies 
{dOnk/dx) = 2TTl/d, where I is an integer. By continu- 
ity of the phase with k, we expect I to have a fixed value 
for a given band, and (vg) = {h/m){k + G) where G is 
some reciprocal lattice vector. For the lowest band, we 
show in Appendix IXI that G = 0. Thus, we arrive at the 
somewhat surprising conclusion that (vs) = Hk/m. We 
suspect that similar results apply in higher dimensions 
but have not been able to show this explicitly. 

The average superfluid velocity should be distinguished 
from the velocity determining the average current density 



1 

Vjv 
n 

nv„k ■ 



^ / — V$„k 



d^rw:^i^{p + hk)wnk 



This velocity is given by [2C 



1 



v„k = TVke(n, k) , 



(9) 



(10) 



where the operator L is deflned as 



L= — 



2m 



Kpt + 25r|4>„k| -A^nk- (14) 



Each distinct solution labelled by the index i corresponds 
to a collective excitation of the condensate and Ei repre- 
sents the excitation energy of the mode. The orthonor- 
mality of the quasiparticle amplitudes is specified by 



X 



d r [u*Uj - v*v.j\ = 5ij . 



(15) 



Since the operator L has the translational symmetry 
of the lattice, the Bogoliubov equations admit solutions 
of the form 



Vi{Y) 



i(q+k)-r 



„i(q-k)-r 



where Mi(r) and Ui(r) have the periodicity of the lattice. 
These functions satisfy 

L^,^-kV,{r) - gnw*^^{v)ui(T) = -EiV,{y) (16) 

with 



where e(n,k) — Etot/N is the energy per particle in the 
state characterized by the mean density n and quasimo- 
mentum k. 

In one dimension, the average current density vanishes 
at the zone boundary k = Tr/d if the mean density is 
below a critical value fie On the other hand, the 

average superfluid velocity is (vg) = hiT/md. This is not 
a contradiction since the local superfluid velocity in Q is 
averaged differently when calculating the average current 
density. 

Dynamical states of the condensate are determined by 
the time-dependent GP equation (Q. For small am- 
plitude excitations, the condensate wave function is ex- 
pressed as 



«'(r,i) = [$„k(r)+^$(r,i)]e- 



(11) 



Lq,k = - — (V-fiq-t-ik)2-fFopt-f25n|w„kr-/Ltnk- (17) 

Our notation emphasizes that q and k play distinct roles 
in the Bogoliubov equations: the former characterizes 
the Bloch-like character of the quasiparticle amplitudes 
while the latter corresponds to the quasimomentum of 
the condensate wave function. 

In the following, we shall also make use of the Hamil- 
tonian 



;ik(q) = - — (V+iq+ik)2-fV;pt-Hgn|w„kp-M«k (18) 

which for q = is just the Hamiltonian determining the 
time- independent condensate wave function ui^k- 
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III. PHONON DISPERSION FOR A 
STATIONARY CONDENSATE 

We begin by considering the simpler situation in which 
there is no superfluid flow (k = 0). In this case, 
the ground state solution of the time-independent GP 
equation can be taken to be real, and "-w^^q j.^g(r) — 
^^n=o k=o('") — '^c(r). It is then convenient to introduce 
the functions '\l}'f = Ui±Vi and to combine the Bogoliubov 
equations into a single equation for |2lj |. 



hlib+ + 2gn,hoiJ+ = Eft. 



(19) 



where /ig is the Hamiltonian 



^o(q) = -7^(V + iq)2 + K,pt+5n|u;o.op-Mo.o, (20) 
zm 

in which the mean field, g7i|wo.oP, of the current-free 
condensate appears. 

To solve (|19|) . we introduce a complete set of Bloch 
states which are solutions of the equation 



^o(q)w„q = e„(q)w 



nq • 



(21) 



This is a linear Schrodinger equation but the solution for 
q = and n = coincides with the self-consistent GP 
solution wofi- By definition, the band energies, e„(q), are 
referred to /J,o,o, so that eo(0) = 0. The functions w„q 
satisfy the periodicity property u'„q(r + R) = i/;„q(r) and 
the orthonormality relation 



(22) 



In addition, at q = they are chosen to be real. 

Since ij^f is itself a Bloch function, it can be expanded 

as 



(23) 



The label i represents a band index m and the Bloch 
wave vector q. However in the following, we will only be 
interested in the lowest excitation band and will there- 
fore drop the label for convenience. Substituting this 
expansion into (|19|l we obtain 

^ A/„„,(q)e„,(q)c„,(q) = £;2(q)c„(q) (24) 

n' 

where 



M„„/(q) = — / d''ru;*q(r)2.gnc(r)u'n'q(r) -I- e„(q)5„„' • 
J 

(25) 

We have displayed explicitly the q-dependence of all the 
variables. 

For a cubic lattice, we anticipate a particular eigen- 
value E{q) which has a linear dispersion of the form 



E{q) = hsq - 



(26) 



Our objective is to derive an explicit expression for the 
Bogoliubov sound speed s. From Eq. l(T^ it is clear that 
in the g — > limit, the eigenvector corresponding to this 
particular eigenvalue will be 



c„(0) = SnO ■ 



(27) 



where n — labels the lowest Bloch band, since only this 
band has a vanishing energy (£o(0) = 0). As a function 
of q, the lowest band energy behaves as 



2mo 



+ Oiq') 



(28) 



which defines the effective mass mg of this band. We em- 
phasize that this band mass is determined by the linear 
Schrodinger equation (|21|l . More will be said about this 
later. The phonon eigenvector is a continuous function 
of q and, as we shall see, behaves as Cn{q) = Sno + C(?^) 
for small q. 

To obtain an expression for s we separate the n — 
equation in (|24ll 

E'^coici) = A/oo(q)eo(q)co(q) 

n'=£0 

from the n ^ equations 

E^Cniq) = M„o(q)eo(q)co(q) 



-E 



^ Mnn' (q)£n' (q)c„' (q) . 



Since eo(q) and E"^ are both proportional to g^, the latter 
equation shows that c„(q) cx q^ for n 7^ 0, as claimed. 
Thus, to order g^, these equations can be replaced by 

i;2co(0) =Moo(0)eo(g)co(0) 

+ E ^^On'(0)£„'(0)Cn'(9) 



= M„o(0)eo(g)co(0)H 
Solving for i?^, we obtain 



.(0)£„.(0)c„,(g) (29) 



E' 



J2 MoniM-^)nn'Mn'0 
nn' 



(30) 



where all quantities within the square brackets are un- 
derstood to be the q = values. The prime on the 
summation indicates that the terms n =^ and n' = 
are excluded from the sum. The matrix M„„' is the ma- 
trix obtained by deleting the first row and first column of 



M„ 



We note that this combination of matrix elements 



can in fact be written as 



(M-i)oo 



(31) 
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Thus we find that the square of the sound speed is given 
by 



Thomas-Fermi Limit 



2mo(M-i)oo 



(32) 



We next relate the sound speed to variations of the 
chemical potential with mean density. Writing for sim- 
plicity wq = wq^o and /io = /^o,o, we have 

- T^V^Wo + VoptWo + gnwl = hoWq . (33) 
The derivative of this equation with respect to n is 

ho{0)wo^n + gwl + 2gnwlwo,n = Mo,nWo (34) 

where we use the notation (• • ■)^n to denote a derivative 
with respect to n. Taking the inner product of H34|l with 
Wq and noting that hoWQ = 0, we find 



A*0,n = Q / ^0 ~?r / '^^^ ^0^0, n ■ (35) 

To solve (|34|) for wo,fij we note that the normalization 
condition in H22|l implies 



(frwo^fiWo = . 



(36) 



Thus, wo,ji is orthogonal to wq and has the expansion 

(37) 



in terms of the (real) q = Bloch functions w„ = 'Wn,q=Q- 
Substituting this expansion in (|34|l yields 



(38) 



Using the expansion (|37|l for wo.n in H35I) with the expan- 
sion coefficients defined by (|38() . we find 

MO,n = ^ ( Moo - Y^'^'^On (m-^) Mn'O ) ■ (39) 
\ nn' / 

Comparing this with H31|l . we see that (|32|) is equivalent 
to 



Too 



n (9^0,0 



Too 



dn 



(40) 



This result for an optical lattice was first given by 
Menotti et al. on the basis of general dynamical con- 
siderations. We see here that it follows directly from the 
Bogoliubov equations and also applies in the case of a 3D 
lattice with cubic symmetry. The small-q expansion can 
be viewed as a systematic way of implementing the slowly 
varying ansatz used by Kramer et al. 13] . The expression 
for s has the same form as for a homogeneous gas, with 
Too replacing the bare mass to and the the density deriva- 
tive of the chemical potential, fJ-o,n, replacing the inter- 
action parameter g. In other words, at long wavelengths 
the condensate behaves as a gas of particles of mass too 
with a compressibility, k, given by = h{d^o^o/dfi). 



The Thomas-Fermi (TF) approximation is valid when 
the density varies in space on a length scale much larger 
than the local coherence length ^ = ^/W/2mgn. In this 
situation, the density is well-approximated by 



"-o(r) ^ -i^J■o- K)pt(r)) 

g 



(41) 



except in regions where Vopt — /^o- This does not occur 
if ^0 > [Vopt] 1 



and we then have 



Ho ^ gn + Vopt 



(42) 



where T^opt is the mean value of the optical potential 
in the unit cell. Thus, /io,n = g as for a homogeneous 
gas. Since the effective potential, V^pt -I- gno, in the GP 
equation is a constant in the TF limit, we would expect 
the band mass, too, to be close to the free particle mass, 
TO. One can in fact show that the deviation of toq from 
TO is proportional to V'q^(C/'^)^j where Vq is the amplitude 
of the potential modulation. Since we are assuming that 
^/d <C 1, Too ~ TO and the TF sound velocity is 



(43) 



stf — \ — 



as for a uniform gas. It should be noted that this result is 
valid even when the amplitude of the density modulation 
is of order the mean density h, provided only that the 
inequality ^/d <^ 1 is everywhere satisfied. 

If fiQ < [K)pt]maxi the Thomas- Fermi density develops 
'holes' in regions where Vopt > Mo- For a one-dimensional 
modulation the density is disjoint, as is the case in two or 
three dimensions for sufficiently small density. In this sit- 
uation long wavelength propagating phonon excitations 
cannot exist within the TF approximation since the nec- 
essary fiuctuations in the number of atoms from one lat- 
tice cell to the next cannot occur. In reality, the GP 
density in regions where Vopt > Mo is small but finite and 
phonon-like excitations continue to exist. However, in- 
creasing the localization of the density in the potential 
minima leads to larger effective masses and eventually 
the sound speed s tends to zero. This behaviour cannot 
be described within the TF approximation. 



B. Weak-Coupling Limit 

The g-dependence of /io,B appears explicitly in (|35|l and 
implicitly through the wave function wq which satisfies 
(|33|l . To extract the dependence in the limit g ^ 0, we 
expand the wave function as wq = w^^ +g{dwQ/dg)g=Q + 
■ ■ •. Since wq depends on g through the combination 
gn, we see that g{dwo/dg) — n{dwo/dn). Thus, the 
combination nwo^n appearing in (|35|l is proportional to 
g in the small-g limit and the second term on the right 
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hand side of H35() is of order . We then have 



Mo,' 



n 

n 



n 



(44) 



As discussed in Ref. [l^, the first term accounts for the 
effect of the lattice on the compressibihty, k, which de- 
creases with increasing localization of the wave function 

Wq"' . The second term shows that /io,n deviates from a 
linear dependence on g as the strength of the interaction 
is increased. An explicit expression for this quadratic 
correction can be obtained from the equivalent expres- 
sion for /io.n in 1)39(1 in which the two terms respectively 
correspond to the two integrals in (|35|l . From the defini- 
tion of the Mnn' matrix in (|25|l . we see that 



2n 



2n 



-M, 



On 



Thus, 



S,vn 



-(0) 



■9 / j3 / (0)\4 



+ 



(45) 



Since the excitation energies are positive, we see that fj,Q^n 
has a negative curvature, which agrees with the numerical 
results in Ref. The interatomic interaction has the 
effect of increasing the width of the wave function which 
counteracts the localizing effect of the lattice potential. 



where the superscript here denotes the order in Vq. The 
properly normalized wave function in the absence of the 

potential is w^^ — 1 and /ip"'' — gn. To first order in Vb, 

/1q^^ — and 



(1) 



-(0) 



-I- 2gn 



■ cos{Gz) 



(47) 



where e 



(0) 



h^G^/2m. In calculating the second or- 

(2) 

der contribution, fi^ , to the chemical potential, the sec- 
need not be determined. 



ond order wave function, 



„(2) 



but the normalization condition /^^^2 '^^ 



„(0),„(2) 



'''0 



— (1/2) j'_^^^dz{wQ ) is required. Thus, the chemical 
potential correct to second order in Vb is found to be 



Mo = gn 



^0 



2(4") 



2gnY 

Taking the derivative with respect to n, we have 



(0) 



Mo,n = g 



-^gy^e^G 



2gnf 



(48) 



(49) 



The weak-coupling limit of this result to lowest order in g 

is A*o,ri = ff(l + 2(Vb/eQ))^). This can be shown to agree 
with the expansion of the first term in (|45|) to second 
order in Vb- 

To complete the calculation of the sound speed, we 
require an equivalent expression for the effective mass. 
This can be obtained by solving 



2m \ dz 



h iqj wog + Vb cos{Gz)woq 
+ {gnwl ~ fJ.o)woq = £o{q)woq , (50) 



where wq is the solution of (|46|) . Since the correction to 
the effective mass is second order in Vb, it is sufficient 
to consider the mean-field potential {gnwQ — /io) to first 
order in Vq. Wc must then solve 



Weak Potential Limit 



fi^ f d 

( ^ + «g ) woq + Vg cos{Gz)woq = eo{q)woq , 



It is also of interest to obtain an expression for the 
sound speed in the case of a weak optical potential 
where perturbation theory applies. For simplicity we con- 
sider a weak one-dimensional periodic potential V^pt = 
Vbcos(G'2;), where G = 2TT/d, applied to an otherwise 
three-dimensional system. The relevant GP equation is 
now one-dimensional. 



2m dz"^ 



Vb cos{Gz)wq 



gnwl 



(46) 



In treating the optical potential as a perturbation, we 
expand the wave function as wq 
and chemical potential as fj,o = 



- ,„(0) J_,„(1)_L,„(2) 



(0) 
Mo 



'^O 



'^o 

(2) 

Mo 



where Vq = Vbe^^'/l^G^+^S"-)- Again treating Vq pertur- 
batively, the effective mass for the lowest band is found 
to be 



1 

Too 



1 

TO 



2Vq^ 



(4^ + '2gny 



(51) 



Inserting (|51|l and H49|) into (|40l) , and discarding the quar- 
tic term in Vb, we arrive at the following expression for 
the sound speed 



\ 



gn 



4gnF(2 



2gnf 



(52) 
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When Eg <^ gn (or ^/d ^ 1 /27r) , this expression reduces 
to 




(53) 



in agreement with the approximation to the sound speed 
obtained by Berg-S0rensen and M0lmer |^. 

To make contact with Ref. 0| , we introduce the recoil 
energy Ej^ — fi^n^ /2mcP — e^^ /i and define 2Vo = (jEr 
(the parameter a is cahed 's' in Ref. IT^). Eq. ^ can 
then be rewritten as 



7(7 



So 



128(1 + 7/2)3 



(54) 



The eigenstates of ho{q) are now labelled by the set 
of quantum numbers {q,n,m,v}, where n is a one- 
dimensional band index, m is the azimuthal quantum 
number associated with the z-component of angular mo- 
mentum, and V labels the different radial excitations. Of 
interest here are axially symmetric solutions (m = 0) 
since these have the character of the phonon mode of 
interest. 

The analysis after ifT^ is followed step by step, the 
only change being the integration volume used in the 
normalization of the states. We thus find that the sound 
speed is given by 



mo d\ 



(58) 



where sq — \J gnjm is the sound speed for the homo- 
geneous gas and 7 is the ratio gn/En — (d/7r^)^. This 
shows that s decreases quadratically with the strength 
of the optical potential, which is consistent with the nu- 
merical results in Ref. This expression is valid if 

2Vb <C or cr <C 4, however, this constitutes a rather 
limited range of the values of a of physical interest. 



D. Radially Confined Condensates 

As a final application of the results derived in this sec- 
tion we consider a condensate that is confined in the ra- 
dial direction. To be specific, we assume a potential of 
the form 



F(r) = Kpt(^) + ^±(p) 



(55) 



where Vj_(p) = raijJi_p'^ /2, that is, harmonic confinement 
in the radial (p) direction. The optical potential is peri- 
odic in the axial direction with periodicity d. This po- 
tential approximates the situation of a long cigar-shaped 
trap with an axial standing light wave. 

Although the geometry is quite different from that con- 
sidered earlier, the previous analysis can be carried over 
with minor modification. The ground state GP solution 
has the property <I>o(p, z d) — ^^{p^z), and has the 
normalization 



1 



d/2 



-d/2 



dz / dVi|$o(r)|2 = A, 



(56) 



where mo is the effective mass of the lowest band (n = 
and V — Q). 

Eq. H58|l is of course valid in the absence of the optical 
potential. Treating the condensate in the TF approxima- 
tion, we find A = ttpq/ gmuj'^ and dpa/dX = gmuj\/2'KiiQ. 
This gives a sound speed 



.9?^o(0) 
2m 



(59) 



where no(0) — po/g is the density at the centre of the 
trap. This result was first obtained in Ref. using a 
different method. In the weak coupling limit, (|45|l ap- 
plies. In this case, the condensate wave function is a 
gaussian. 



(0) 



exp 



2h 



and one again obtains the result in (|59|l for the sound 
speed, as found previously [2^. 

When the optical potential is strong the condensate 
becomes localized on each site. In this situation, the 
tight-binding approximation is a useful method for deal- 
ing with the system 10. 11. 24j . Approximating the con- 
densate wave function as $(r) = X]iCi/i(r) where /^(r) 
is a function localized on the i-th site and normalized to 
unity, the energy of the system is given approximately by 



Eu 



]£o\ci\'' 



t^iCiCi+i + c,c: 



which defines the mean linear density A along the length 
of the condensate. As before, it is convenient to define 

$o(i-) = VXwo{r). 

The Bogoliubov excitations in the present situation 
have a Bloch wave character along the axis and are ob- 
tained from (|19|l with the Hamiltonian 



d 

oz 



+ V + gX\wo\^ - fio ■ 
(57) 



where Eq is an on-site energy, t is a hopping matrix el- 
ement connecting the amplitudes on nearest-neighbour 
sites and g is an effective interaction strength. For the 



where ^ is the number of atoms 



ground state, |cj| : 
per site. We then have 



-Btot = (eo - 2t)N + -guN . 

Assuming that the parameters Eq, t and g are den- 
sity independent in the extreme tight-binding limit, we 
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have ^0 = dEtot/dN = {sq — 2t) + gv and dfio/dX = 
gd. Within the same approximation the band energy is 
e{q) = eo + 9>^ — 2t cos{qd), which gives an effective mass 
ma — fi^ /2td^. Thus the BogoHubov sound speed from 
()58|l is stb — "Igvtd^ / 'h? , which is the result obtained 
by Javanainen [lc| . 



IV. PHONON DISPERSION FOR A MOVING 
CONDENSATE 

We turn next to the derivation of the phonon disper- 
sion relation for the case where the condensate is "flow- 
ing" through the optical lattice. We address this problem 
by directly solving the Bogoliubov equations in 1)16(1 for 
a specific condensate wave function Wn\^ in the limit of 
small wave vectors q. The structure of these equations 
is quite different when k ^ and the method used in 
the previous section to determine the dispersion relation 
can no longer be applied. In fact, the analysis is much 
more intricate as will soon become apparent. The re- 
sults we obtain confirm the more intuitive hydrodynamic 
approach presented recently which describes the dy- 
namics of the system in terms of slowly varying hydro- 
dynamic variables (density and momentum). By includ- 
ing small length scale variations, our approach in a sense 
provides a "microscopic" derivation of the hydrodynamic 
equations that one would expect to be valid in the long 
wavelength limit. 

Since the solution in the small-q limit is required, we 
rewrite so as to display the q-dependent terms ex- 
plicitly: 



&n(q)Wn,-k 



— • (p + Tik) + —— 
m Zm 



,,2 



^•(p-^k) 

m 



—gnwQj^v = Eu 



2m 



where 



(61) 



(62) 



According to this definition, egk = 0. The functions 
are an orthonormal set with normalization given by H22|l . 
Although we use the same notation, it should be noted 
that these functions are distinct from those defined in 
l|?T|l . In addition, we may assume w„,_k(i") = w*^f^{r) 
and e„,_k = £„k- 

Substituting these expansions into (|60|l . we obtain the 
matrix equations 

£nk + -7: On + / — q • -r rm'fln' 

2m J '—f m 



+ (A„„/a„' - Bnn'bn') = Ea„, (63) 



£nk + -Z I On 

2m 



+ J2 iKn'K' - Bl^,a^,) = -Eh,,, (64) 



where we have defined the matrices 

Ann'O^) = 7^ I 'w'^kgn\wQu\'^Wn'-kd^r, (65) 

-Bnn'(k) ^ <k5^W0k<'kC^^''> (66) 



p , 



-gnwolu 



-Ev, (60) 



(k) = i^<k(P + ^k)t.„,kdV (67) 

Due to the inversion symmetry of the lattice, the Bloch 
states at the zone centre (k — 0) can be chosen to be 
simultaneous eigenstates of parity. Although parity is 
not a good quantum number for states with nonzero k, 
each band can nevertheless be assigned a parity index 
r]n = ±1 such that 



where p — {h/i)V is the momentum operator. The GP 
Hamiltonian here is /ik(q = 0) as defined in (|18|l . In 
these equations, we have adopted the index n = for 
the condensate wave function. We will usually think of 
this state as the lowest Bloch state solution of the GP 
equation, although it in principle could correspond to an 
arbitrary excited band. For simplicity we have dropped 
the index on the quasiparticle amplitudes u and v and 
the excitation energy E as we will only be considering 
the phonon-like excitation. 

It is clear that for q = H60() admits a solution with 
u oc wok, ^ oc Wq[^ and E ~ 0. For finite q we seek 
solutions in the form of an expansion in eigcnfunctions 
of /i±k, namely. 



E 



w„ _k(-r) ?7„w„k(r) , 



(68) 



This property is proved explicitly in one dimension in 
Appendix ^ and can also be shown to follow to lowest 
order in k by means of k-p perturbation theory. Together 
with the conjugation (time-reversal) property w^^{r) = 
Wn.-kir), we have 



(69) 



This important property is used throughout the following 
discussion. For example, it can be used to show that the 
matrices in (|65I67|I satisfy the following relations 

A;„,(k) = A„.„(k) = A„„,(-k) = ^nr]„'An„'0^) 

S*„,(k) = B„„/(-k) = rjnTjn'Bnn'i^) = ?7„?7„' B„/„ (k) 



Prm'(k) — Pn'n(k) 



Pfin' ( k) — P„„' (k) 



(70) 
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Note that A and P are hermitian while B is not. In 
addition, we see that A„„/ (0) and B^n' (0) are real and 
nonzero only for pairs of states having the same parity 
index, while Pnn' (0) is purely imaginary and only couples 
states with opposite parity. 

In solving (|63|l and (|64|l . it is convenient to define the 



following linear combinations: 



i(a„ + 77„6„) and 



dn = hio-n — Vnbn)- Introducing these variables into H63() 
and Hd4I) . and making use of the relations in H7U|) . we 
obtain the equations 



V ^ ~2m ) ^ 2^ [^Ann' Bnn 
^ ^ n' 

Eh 
— q • P„„'d„' Edn, 
m 



(71) 



Eh 
— q • Pjm'Cn' = ECn, 
m 



(72) 



where we have defined the hermitian matrix Bnn' ~ 
Bnn'Vn'- As in our earlier analysis, we anticipate that 
E{q) will depend linearly on q for q ^ 0, but due to 
the quasimomentum k of the condensate, it no longer 
depends simply on the magnitude q. To extract this de- 
pendence we systematically expand the coefficients c„ (q) 
and (i„(q) as a series in powers of q. Specifically, we write 



c„(q) ^ Mq)(cl°' + c« + c(^)---) 
rf„(q) = Mq)(4°^+4'^ + 4''---) 



(73) 



where the superscript indicates the order of q in the 
respective terms (here, order signifies similar powers of 
the vector magnitude q H^). The factor h{q) contains 
the nonanalytic behaviour of the coefficients which is re- 
quired in order to satisfy the normalization condition in 
(|15|l . For a homogeneous system, h{q) cx q^^^^, and we 
expect a similar dependence in the case of a lattice. In 
the following, it is sufficient to note that this factor is the 
same for both coefficients and can therefore be ignored 
in developing a systematic g-expansion. 

We noted earlier that u cx wqu and v cx u'o,-k for q — > 
which, according to (|61|) . implies that a„(q) oc Sno and 
6„(q) cx Sno in this limit. If the state wok is an even 
parity state (770 = +1)? as we assume in the following, 
we must then have Cn^ — Sno and dn^ = 0. With this 
information, H71|) gives to 0(9°) the equation 



EOk + (AnO — B„ 







(74) 



which is satisfied since £ok — and i3„o — A^o- To 
©(gi), IZH) gives 

£nkci^^ + J2 iAiW - Bnn' )c^n' = « • (75) 



Similarly, l(7^ gives 

V Nnn'd^n' = ^SnO " -q ' FnO , (76) 
n' 

where we have defined the hermitian matrix 



Nnn' (k) = Ann' (k) + S„„/ (k) + £„k(5„ 



(77) 



Eq. H75|l is homogeneous and indicates that ci^'' = 0. On 
the other hand, H7t)|l can be solved for dn^ in terms of the 
unknown excitation energy E. To determine the latter, 
we must consider (|71ll to 0{q^), obtaining 



2„2 



ha ~ 

Enkci^' + -:^SnO + X1(A„„- - Bnn')c''n, 



(2) 



2m 



E 



-q ■ P nn' diV ^Edi'K (78) 



m 



Setting n = in this equation, and noting that eok = 
and that Aon — Bon, we find 



2m 



(79) 



Since dn"^ is itself linear in E according to (|76|l . we see 
that (|79|l is implicitly a quadratic equation for E which 
can be solved to determine the excitation energy to low- 
est order in q. However, to do so directly would not 
reveal the interesting dependences on various physical 
parameters that in fact emerge. As seen in the k = 
analysis, the excitation energy could be related to the 
variation of the chemical potential with mean density n. 
This remains a quantity of interest in the present case, 
but we must also consider variations of the chemical po- 
tential with k. We thus turn next to the determination 
of A*ok,fi = dfiok/dn and ^ok,i = 9^ok/<9fci. 



A. Determination of ^ok.n 

The chemical potential /iok is determined by the self- 
consistent solution of the GP equation, /ik^ok = 0, with 
hk given by H18|) with q = 0. The variation of this equa- 
tion with respect to h gives 

hkWok^n + g |wok|^ Wok + gn (wok^^Wok + WokWok.fi) Wok 

= MOk.nWok, (80) 

where wok.fi = dwok/dn. Taking the derivative of the 
normalization H22|) with respect to n, and using (|69|l . we 
find that 

/ WgkWok.ftrf^r = 0. (81) 
Thus, wok.fi is orthogonal to wok and has the expansion 

WOk,n = anWnk , (82) 
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where as before, the prime on the summation indicates 
that the n — term is excluded from the sum. Inserting 
this expansion into (|8()|l . and taking the inner product 
with respect to Wnk, we obtain 

1 .r-^' 

TT^NnO + y ^ Nnn'Ctn' - M0k,n<5„0 = , (83) 

In ^ — ' 



where we have noted that a* = ry„a„ as a result of the 
symmetry property (|68|l . Setting n = in (|83|l . wc find 



1 

/^Ok,fi = 77^ ^^00 + > Nona„ 
Ml. 



2n 

n 

The set of equations for n 7^ has the solution 

1 NT^' 



(84) 



(85) 



where N is the reduced matrix obtained by deleting the 
first row and first column from N. Thus we find that 



^(^00- 5l'^0„(A^~')n»'^n'0 J ■ (86) 
\ nn' J 



This is analogous to and reduces to it in the k = 
limit since the M and N matrices are then the same (note 
that Wno are defined to be real). 



B. Determination of /iok,i 

We follow a similar method to obtain //ok,z- Taking the 
derivative of /ikti'o.k = with respect to fcj, we have 



r h 

m 



MOk,i 



wok 



+/ikWok,i = 0, (87) 



where wok.i — dwQk/dki. Noting the orthogonality of 
each vector component wok.i with wok, we have the fol- 
lowing expansion 



(88) 



where the expansion coefficients, /5j„, define the Carte- 
sian components of a vector /3„. Inserting H88|l into (|87|l . 
and taking the inner product with Wnk, we obtain 

h V — a' 

— (-Pi)nO + > Nnn'Pin' - AiOk,i<5nO = 0. (89) 

m ^ 

n' 

where, as before, we have used = r?„/3m. An expres- 
sion for /iok,i can be found by setting rt = in H89(l : 



MOk.i — — (^j)oO + NQn'Pin' 



(90) 



The set of equations for n 7^ yield the solution vector, 

ft NT-', 



A, 



\ J- ^ 



(91) 



and thus, 
MOk,i — 



^ ((^Ooo-^'iVo„(iV-i)„„,(POn'o] . (92) 

\ nn' / 



C. Excitation Energy 

These results will now be used to obtain an expression 
for the excitation energy E from (|76|l and H79(l . Setting 
n = in (|76|) we have 

V Vo„-4V = £; - -q • Poo - iVoo4'^ (93) 



while for n 7^ we see that 



di'^ = -EV"')nn' f-q • Pn'O + ^n'O^'^) (94) 



The quantities on the right hand side of this equation are 
in fact related to the expansion coefficients a„ in H85(l and 
Pin in (|91|l . We find the simple relation 

=2n4'^a„-Kq-/3n. (95) 
Using this result in H93|l . we have 

(^Nqo + 2nE'iVonan^ = - • Poo 

- EVonq-/3n (96) 



(97) 



which with (|84|) and 190|) can be written as 

2n^ok,nrfo^^ = E -q - VkMok • 

We see that d^^ and E are now related to each other 
through physically meaningful and calculable parame- 
ters. 

We now substitute into lf7S|) to obtain 

ft2g2 , ft 



2m ^-^ m 



J2 -(q-Pon)(q-/3n) 



^ (e--ci-Poo- V'-q • Pon2na„ ) ■ (98) 

\ TO ^ TO / 

With (|85|) . the sum on the right hand side becomes 
V'-Po„2na„ ^ -V'-Pon(^"^)„n'A^„'o 



= -Y'NoniN-Xn'-Pn' 

^ — ^ m 

nn' 



(99) 
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In going from the first to the second hne, we have used the 
fact that all the matrices have the transposition property 
Mn'n = TlnVn' Mnn' ■ With the cxpression for /zok.i in JHOIj 
(|98|l thus becomes 



2„2 



2m 



J2 -(q-Pon)(q-/3„) 
^ — ^ m 

71 



(100) 



Eliminating dg^"* from H97|l and (|l()()|l . we finally obtain 
E = q- VkMok 



A 



Ok, ft 



2m 



+ -y"'(q-Pon)(q- Ai 



.(101) 



The sign of the square root is chosen to be positive to 
give a positive excitation energy in the k — > limit. The 
final quantity to interpret is the summation within the 
square root. 



D. Effective Mass Tensor 



The square bracket in IjlOll) involves the tensor 



1 



1 

m I m 
1 

m 



2 



E'(^')o„(/3,). 



^Mo - All'(^Oon(A^~')™'(Pj)«'o. (102) 



This expression is similar to the usual effective mass ten- 
sor defined on the basis of k • p perturbation theory, al- 
though the structure of the summation is different. To 
make contact with the k • p expression we consider the 
ISSnn' matrix in the k ^ limit. Quite generally, this 
matrix has the block structure 



N(k) ^ ( 



A+ B+ 
A B +D 



(103) 

where the blocks are defined according to the parity in- 
dex of the various states. For example, the block in 
the upper-left corner contains matrix elements between 
states with a positive parity index, = -t-1. The diagonal 
matrix D contains the energy eigenvalues e„k on its diag- 
onal. In the limit k — s- 0, we have B(k = 0) = A(k = 0) 
and Ann' = if r\n 7^ '7n' • Thus, in this limit we have 



N(k = 0) 



2A4 



D, 








D 



(104) 



that is, N(k — 0) is block-diagonal, which of course is 
also true of its inverse. Since, Pi only connects states of 
opposite parity in the k = limit, we thus see that 



lim ( — 

k^o V rn 



1 



^' {Pi)0n{Pj)n0 



(105) 



This is precisely the effective mass tensor obtained by 
means of k • p perturbation theory [23 as applied to the 
Hamiltonian /io(q) in H20|) . The tensor defined in H102|) is 
a generalized effective mass tensor in that it depends on 
the presence of a superfluid flow (k 7^ 0). Also because of 
this, it is no longer diagonal despite the cubic symmetry 
of the optical lattice. 

To complete the identification of (to^^)^, we consider 
variations of the GP equation with respect to the conden- 
sate wave vector k. The second derivative of h^wok = 
yields the equation 

hk,ij'Wou + /ik,iWok,j + hKjWok,i + hkWoic,ij = . (106) 
Here, 

hk,i = —iPi + fiki) + gnc,i - Mok.j , (107) 
m 

hw,i3 = — ^ij + gn^ij - Aiok.ij ■ (108) 

TO 

The inner product of H10t)|) with wok gives 



-f; / ^ok ( ^k.yWok + hk^iWokj + /ikjWok,i d r = 



(109) 



The first integral is 



(110) 

while the integral of the next two terms gives the result 

77 / 1i'0k(^k,jW0kj + /ikjWok,i)c?^?' 



VL 



-E'l 

m ^ — ^ 



{Pi)Qn{Pj)n + (-Pj )on (A )n] 



-I- 



(111) 



We thus find 
-Si 



A'Ok/H — Oij -I > [Pi)On{Pj)n 

m m ^-^ 



2nn 



(112) 

With this result we see that the effective mass tensor 
defined in l|102|l can be expressed as 
1 92 



,■■ h? dkidkj 



2nn 



nid-^r] . (113) 



E. Relation to Energy Density 

This last result can be related to the total energy Q 
in the state <i>ok- Defining the energy per particle as 
Etot = Ne{n, k), we have 

?(n,k) = ^ / w*^(~^{V + ik)^ + Vop?\wokd^r 



2m 



(114) 
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Comparing this with /iQk, we see that 



e(n, k) = /iok ^ 1^ ^ Iwowl^d^r 



(115) 



which is the expression in brackets in pi3|) . Thus the 
effective mass tensor is given by 



1 d^e 
dkidkj 



(116) 



For a cubic lattice, /^ok has an expansion of the form 
Mok = /^oo + h^k'^ /2ra^ + Similarly, e(n, k) = 

e(n, 0)+/l^/c^/2me + - • •. However, as proved by pU5|) . the 
parameter mg is in fact the band mass mg defined by the 
Hamiltonian (|20|l . In other words, the correct effective 
mass parameter can be extracted without solving the GP 
equation for wok (with k ^ 0) self-consistently. We note 
in passing that direct differentiation of (|115|l establishes 
the relation (js) = '^k{ne)/h |20|. 

With these results, the phonon energy H1U1|) can be 
written in a compact form. Defining the mean energy 
density as e = ne, the Bogoliubov excitation energy at 
long wavelengths is given by 



(117) 



where we use a repeated summation convention on the 
Cartesian indices i and j. This is precisely the expression 
given by Machholm et al. who argued that the dynam- 
ics of the system at long wavelengths could be based on 
a hydrodynamic analysis. Since their approach arrives at 
(|117|l in a more economical fashion, it is useful to sum- 
marize the essential assumptions on which it is based. 

The central assumption is the existence of an average 
phase fluctuation, {6{r,t)), that varies slowly in space 
and time. Expanding this average phase as 

(0(r + Ar,t + At)) = (0(r,t)) + V(0(r,i))-Ar 



dt 



-At + ---, (118) 



one identifies V(0) with the local wave vector, (k), and 
d{9)/dt with —{ij)/h where (/i) is the local chemical po- 
tential. The equation of motion for the local wave vector 
is thus 



(119) 



The second hydrodynamic equation is the continuity 
equation 



dt 



V • (j,) = , 



(120) 



where (js) is the local current density. The current den- 
sity and chemical potential are then assumed to be given 
by the usual expressions for a uniform system, namely. 



1 , ^ 

TV(k)e, (a*) 



de 
d{n) 



(121) 



where e((n,), (k)) is the average energy density for a uni- 
form optical lattice, viewed as a function of the local 
density (n) and wave vector (k) By expanding the 

variables as {n) — n-\-5n and (k) = k-|-(5k, one obtains a 
pair of equations for the fluctuations which admits wave- 
like solutions with frequency lo = E/h and wave vector 
q. The dispersion relation found is identical to H117|l . 

It is clear that the assumptions made in the hydro- 
dynamic approach are completely justified. The average 
energy density e is the fundamental quantity determining 
the excitation energy at long wavelengths, as confirmed 
by our systematic g-expansion. The additional informa- 
tion provided by the expansion technique are the pertur- 
bative expressions for 9/xok/9n, Vk/^ok and {l/m)ij as 
given by ((HSl, (EU and H102|) . respectively. 



F. Discussion 

For small k, e{n, k) ~ e{n, 0) + nh^k"^ /2mQ 
the Bogoliubov excitation energy is 



E/h ~ ftq • k 



d_ 

dn 



n 

Too 



sq 



and 



(122) 



where s is the sound speed for the condensate at rest. 
This result was given previously by Kramer et al. 1 1 8| . 
The energy first becomes negative when the superfiuid 
flow satisfies k > TO^s/fi, where — d{n/mo)/dn, 

which defines the Landau criterion for energetic instabil- 
ity at long wavelengths in an optical lattice. The region 
of energetic instability was mapped out for arbitrary q by 
Wu and Niu ^,2^ and Machholm et al. JJ. In this re- 
gion the energy of the superfiuid state is no longer a local 
maximum. As a result, transitions to lower energy states 
can occur spontaneously provided a means of conserving 
energy and quasimomentum is available. 

The excitation energy given by pi 7(1 becomes imagi- 
nary when the argument of the square root is negative. 
This signals a dynamic instability whereby the amplitude 
of the condensate fiuctuation grows (or decays) in time. 
Of the two factors in the square root, e^ij, or equivalently 
the effective mass tensor, (I/to),^, is the most physically 
relevant. 

It is instructive to examine the latter in the weak 
potential limit. We consider for simplicity the one- 
dimensional situation discussed in Sec. IIII CI Repeating 
the perturbative analysis in Sec. IIII CI for the case k 7^ 0, 
we find 



e{n,k) = ^gn- 



^ I. 



^0 



2 ( -I- 2gn - As 



(0) 



(123) 



We note that this expression becomes singular at a wave 



vector kc satisfying e^^ -I- 2gn — Ae)^'' — 0, which gives 



(0) 



fee 



(124) 
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where fco = mso/h. The singularity is indicating the 
breakdown of nondegenerate perturbation theory, but 
provided k is not too close to kc, we can use (jl23ll to 
evaluate the effective mass. To second order in Vq we 
have 



Too (A:) 



1 

TO 



1 - 



4°) + 2gn - 4e^' 



(0) 



(0) 
G 



2gn 



(125) 



At k = 

For TO, 
Afc = > 



we recover the result (|51|l found in Sec. IIII CI 
^ to go to zero, k must be close to kc- With 
, — k, we find 



Afc 



Vn 



2gn 



2/3 



(126) 



that is, the wave vector fc approaches fcc as Vq — > 0. We 
note that at fc = fcc ~ the perturbative correction to 
the energy in H123|) is still small so that the perturbation 
theory estimate of where TOq ^ goes to zero is reasonable. 
We thus expect a dynamic instability to set in when fc ~ 
kc in the weak potential limit. 

This condition for the dynamical instability is the 
q = limit of the result given in Refs. 0| and [7|| . The Bo- 
goliubov excitations of wave vector q in an homogeneous 
gas with current (js) — nhk/m have the energies 



E±{q) 



± 



h^kq ngh^q^ h'^q'^ 



Am? 



(127) 



We follow Wu and Niu 6] in referring to the modes with 
the plus (minus) sign as phonons (anti-phonons) . The 
former correspond to physical excitations in that their 
normalization is given by l|15|l . The effect of an optical 
potential is to couple an anti-phonon mode with wave 
vector g to a phonon mode with wave vector q — G. The 
condition that E- (q) = E+ {q — G) implies that the two 
modes are resonantly coupled and gives the critical wave 
vector 



fc§(G 



(128) 

For q — Q, this gives the critical wave vector in H124|) . The 
expression in (|128|l was shown in to account for the 
boundary of the dynamically unstable region in the weak 
potential limit. In fact, it can be shown by means of de- 
generate perturbation theory f Appendix IB|| that impos- 
ing a weak optical potential indeed gives rise to complex 
Bogoliubov eigenvalues. 

Alternatively, the condition E-{q) — E+{q — G) can be 
written as E+{—q)+E^{q—G) = 0. This was interpreted 
by Machholm et al. as a Landau criterion for the emis- 
sion of two phonon excitations with zero total energy. 



Although this physical interpretation is appealing, it is 
not clear how it can be used to actually determine the 
rate at which the excitations are being produced, short 
of performing the perturbation analysis carried out in 
Appendix IbI in terms of phonon and anti-phonon modes. 

We thus see that the phonon-anti-phonon resonance 
condition, or alternatively the two-phonon Landau crite- 
rion, is consistent with the effective mass condition for a 
dynamical instability in the q ^ limit. A similar state- 
ment can be made in the weak coupling limit {g 0). 
Wu and Niu @ noted from their numerical analysis that 
one boundary of the dynamically unstable region is given 
by the condition eo{q + k) — eo(fc) = £o{k) — eo(fc — q) 
where eo{k) is the band energy for the optical potential 
by itself. In the small-g limit, this condition becomes 



afc2 



0. 



Thus the onset of dynamical instability at g = in the 
weak coupling limit is again given by the point at which 
the inverse effective mass goes to zero. 



V. CONCLUSIONS 

We have studied the long wavelength phonon excita- 
tions in a three dimensional optical lattice. By making 
use of a systematic expansion of the Bogoliubov equa- 
tions in terms of the phonon wave vector q, we obtain 
the phonon dispersion in the long wavelength limit. Our 
result (|4()|l for the current-free state defines the sound 
speed in terms of the effective mass toq and variations of 
the chemical potential with n and agrees with the result 
given by Menotti et al. 17]. The effective mass is de- 
fined quite generally in terms of the energy per particle, 
e(n, k), but can also be calculated using the current-free 
GP Hamiltonian in the k limit. We present analytic 
expressions for the sound speed in the Thomas-Fermi, 
weak potential, weak coupling and tight-binding limits. 

For the current-carrying case, we rederive the disper- 
sion relation obtained by means of a hydrodynamic anal- 
ysis (see also 0). Our approach confirms that the 
dynamics at long wavelengths is defined by the local en- 
ergy density e((n), (k)) viewed as a function of the slowly 
varying local density, {n{r)), and local condensate wave 
vector, (k(r)). At long wavelengths, dynamical insta- 
bilities arise at the point where the generalized effective 
mass tensor has a vanishing eigenvalue. 



APPENDIX A: REFLECTION SYMMETRY OF 
BLOCK FUNCTIONS 

In this Appendix we give a proof of the symmetry prop- 
erty (|68|l used throughout our analysis. We do this for 
the one-dimensional case for which the wave function is 
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a solution of 



2m da;^ 



(Al) 



where the potential is periodic, V{x + d) = V{x), and is 
assumed to have inversion symmetry, V{—x) — V{x). In 
the context of the GP equation, V{x) = Vopt{x) + gndx) 
and the inversion property is valid if the condensate den- 
sity also satisfies nc{—x) — ndx). This is ensured if the 
wave function has the property we wish to prove. 

We seek solutions of the Bloch form, il}{x + d) = 
gifcd^^^;). Due to the inversion symmetry, the linearly 
independent solutions of ljAl|l can be chosen to be even 
[ipe) or odd (-00 ) functions of x and '4'{x) can be expressed 
as the linear combination 

ip{x) = ail}e{x) +b-iljo{x) . (A2) 

The periodic part of the Bloch function is then 

Wk{x) = e-*^-" {aij,{x) + bM^)) ■ (A3) 

The two independent solutions at energy E are chosen to 
have the normalization 



1 



d/2 



d/2 



\lJje,o{x)\'^ dx = 1 . 



(A4) 



Imposing the Bloch condition, we obtain the relation 



^ = -^tan^r^ 

iPe Wo 



(A5) 



where all the functions are evaluated at x = d/2. Since 
Tpf. and Tpo are functions of the energy, E, this equation 
determines the band energy Ek- Clearly E^k ^ Ek- If 
Eq is the band energy at the zone centre (fc = 0), we 
must have either ip'^[d / 2; Eo) = or tpo{d/2;Eo) = 0. 
The former defines what we shall refer to as an even- 
parity band, while the latter defines an odd-parity band. 
The small-A; behaviour of E is thus readily obtained from 
these properties. For example, for an odd-parity band we 
have 

Md/'2;E)^ME-Eo) + --- (A6) 
where ipo = d^lJo{d/2] E) / dE\E=Ea- We then have 



Ek = Eq 



(A7) 



The coefficient of fc^ defines the effective mass of the 
band. A similar result applies in the case of the even- 
parity bands. 

Once the energy eigenvalue for a given k is known, the 
coefficients a and b are related by 



b Tpe f kd 
— ~ i — tan — 
a Tpo \ 2 



(A8) 



For a given band, n, the ratio b/a is a continuous function 
of fc. At A: = we choose Wk=oix) to be real and assume 
that it is a parity eigenstate. In this situation, we must 
have either b{k = 0) = (even-parity bands) or a{k ~ 
0) = (odd-parity bands). 

The normalization of Wk leads to the expressions 



1 



1 + A2 ' 



\b\' 



1 + A2 



where 



Xik) = — tan I ^ 



(A9) 



(AlO) 



For an even-parity band A — > as fc — )■ 0, so that a(fc) 
1 and b{k) 0. In this case. 



Wkix) 



VVTTa^ VTTa^ 



(All) 



Since A(— fc) = — A(fc), we see that w^k{^x) = Wk{x). 
On the other hand, for an odd-parity band, A(fc) — > oo 
as fc ^ 0, and b{k) 1. As a result, we have 

^,(,) = e-^- ( + ^^"i , (A12) 



x/TTa^ 



which implies w-k{—x) — —Wk{x). We have thus shown 
that the Bloch functions have the property 



w-k{-x) = ±Wk{x) 



(A13) 



where the positive (negative) sign corresponds to the even 
(odd) parity bands. Together with the conjugation prop- 
erty w1{x) = w-kix), we have w'^{—x) — ±Wk{x). 

For an even-parity band, the real and imaginary parts 
of Wk are 



^Wk{x) 



1 



{cos{kx)'ipe{x) + Xsin{kx)^po{x)) 



3?i'fc(x) = {— sin{kx)'ipe{x) + Xcos{kx)ipo{x)) . 

V 1 + A^ 

Thus the real part is an even function of x with the 
property d?iiwk/dx\±d/2 — 0, while the imaginary part 
is odd and Qwk{id/2) = 0. The opposite is true of 
an odd-parity band. One can show for an arbitrary fc 
in the lowest band that there is no net change in the 
phase Ok{x) = teLii^^{'^Wk{x)/^Wk{x)) as x varies be- 
tween —d/2 and d/2. We make use of this result in Sec.lTTl 
The method described above cannot be used in three 
dimensions, but perturbation theory allows one to infer 
the same symmetry property. We write the Schrodinger 
equation for the Bloch function Wk(r) as 



{ho + SV)w^{r) 



(A14) 



where hn = -{h^ /2m)^'^ + h'^k^/2m + V{r) and SV = 
{Kk/m) ■ V. The eigenfunctions of ho, u)„(r), with 
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eigenenergies are chosen to be parity eigenstates. 
The state Wk to first order in 6V is 



• {n/m)k ■ Pn'n 
1 

En ~ En' 



(A15) 



where P„'„ is the k = momentum matrix element 
defined in ((TDJ. Since this matrix element only cou- 
ples states with opposite parity, we see that w-k(^r) = 
±Wk(r), depending on the parity of the state n. Thus, 
the Bloch state exhibits the symmetry property to low- 
est order in k. It is evident that this argument can be 
extended to higher orders in perturbation theory. 



APPENDIX B: DYNAMIC INSTABILITY IN THE 
WEAK POTENTIAL LIMIT 

As pointed out by Wu and Niu ^ , the boundary of the 
dynamically unstable region in the weak potential limit is 
given by the condition E_{q) = E^{q — G) where E±{q) 
is given by (|127|l . These energies are the eigenvalues of 
the Bogoliubov equations 



with 



Bn = 



ft2 d2 



u± 
v± 



u± 



(Bl) 



(0) 

-2ikz 



-gne 



2ikz 



The phonon mode of interest is 



(0) 



u+ 



i(q-G+k)z 



with 



{eq-G + E+)/2E+ and b+ 



{iq-G - E+)/2E+ where £g_ 



q-G 



JO) 
--q-G 



gn and 

i?+ = E+{q — G) — h^k{q — G)/m. The normalization 
of the mode is — 6^ = 1. The anti- phonon mode 
which is coupled to the phonon mode by a weak optical 
potential, Kpt = Vbcos(G'z), is 



^_^z{q+k}z 
l)_gi{q-k)z 



(B4) 



with a_ 



{Sq - E-)/2E^ and 6_ 



ieq + E-)/2E^, where E- = E+{q) - h^kq/m. 



We note that in this case the mode has normalization 
- = -1. 

The degeneracy of the phonon and anti-phonon modes 
{E^{q) — E+{q — G) = Eq) suggests that we seek a 
solution of the Bogoliubov equations 



B 



in the form 



= A 



B 



(B5) 



(B6) 



Expanding the operator B to first order in the optical 
potential, we have B = Bq + Bi with 



R _ / Vopt + 2gn{wi + wl) -2gne^'''^wi 

^1 ^ I n„-,„~2ikx„,,* 



-2gne 



Vopt + 2gn{wi -f wl) 



Here we have written the condensate wave function as 
^k{x) = -\/He''^^(l + Wl + ■ ■ ■) where the first order cor- 
rection is wi{x) — a+e^'^^ + a_e~ 

-(0) 



a± 



{sff + 2gn. 



4.(0) (0) 

*^G ^k 



(B8) 



Taking the inner product of IjBSp with {u*^_ u^) and 
(ul wl), and noting the different normalizations of the 
two modes, we obtain the matrix equation 



Eo-E A 
A -Eq + E 



0, 



(B9) 



where the real coupling parameter A is given by 



— + 2gn{a+ + a_) ) (a+a_ + 



(BIO) 



A nontrivial solution to the matrix equation is obtained 
if 



E = £^0 ±i|A| 



(Bll) 



Thus, the line in the k-q plane defined by H128|l lies within 
the region of dynamical instability when Vq is finite. As 
emphasized by Wu and Niu , the dynamical instability 
in the weak potential limit arises from a resonant cou- 
pling between phonon and anti-phonon modes. 
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